This RMarkdown explores the comparability of several data products associated with the NPS water balance model (WBM):
We will be using Petersburg National Battlefield (PETE) as our location to compare each of these products’ outputs.
For this comparison, we used the R scripts associated with the WBM found on AWS (https://doimspp.sharepoint.com/:f:/s/nps-nrss-imdiv/EkhiOXJpNKNMpTYuxNO5mBoBqZBScDcbpK5M03Ha61PoCw). These folders contain code to calculate WBM variables for PETE. In the example provided, 10 random points within the park centroid’s MACA grid cell are selected, then used to run the WBM. At the end of the workflow, the 10 separate runs are averaged for each day.
First, let’s load in the functions associated with running the WBM in R:
source("aws_r_tutorial/wbm_functions/ET_functions.R")
source("aws_r_tutorial/wbm_functions/WB_functions.R")
Next, we can load in the metadata associated with the 10 randomly-selected points. This metadata is used as inputs in the WBM and is generated in a separate R script ahead of calculating the WBM.
wb_sites <- read.csv("aws_r_tutorial/PETE_site_parameters.csv")
We also load in the historical GridMET variables associated with our centroid points in PETE. Because we are currently unable to locate the input dataset used in the tutorial hosted on AWS (https://doimspp.sharepoint.com/:f:/s/nps-nrss-imdiv/EkhiOXJpNKNMpTYuxNO5mBoBqZBScDcbpK5M03Ha61PoCw), we are instead using the climate data hosted on AWS for the PETE’s centroid (https://parkfutures.s3.us-west-2.amazonaws.com/maca-tprh-data/index.html). Moreover, we are choosing to start executing the WBM in 1979; we are under the impression that the year 1979 is used as a “warm up” period for both the centroid and the gridded data sets.
Historical <- read.csv("aws_data_downloads/centroid_data/PETE_historical.csv")
Now, we are ready to run the WBM! The code below was pulled directly from the tutorial folders hosted on AWS (https://doimspp.sharepoint.com/:f:/s/nps-nrss-imdiv/EkhiOXJpNKNMpTYuxNO5mBoBqZBScDcbpK5M03Ha61PoCw). In it, the Oudin method is used for calculating daily PET.
n <- nrow(wb_sites)
#Threshold temperature (deg C) for growing degree-days calculation
T.Base = 0
#Method for PET calculation
Method = "Oudin" #Oudin is default method for daily PRISM and MACA data (containing only Tmax, Tmin, and Date).
#Date format
DateFormat = "%m/%d/%Y"
############################################################ CREATE CLIMATE INPUTS #############################################################
#### Historical
# Convert pr.In to mm and F to C
Historical$ppt_mm <- (Historical$Precip..in.*25.4)
Historical$tmax_C <- 5/9*(Historical$Tmax..F. - 32)
Historical$tmin_C <- 5/9*(Historical$Tmin..F. - 32)
Historical$tmean_C <- 5/9*(Historical$Tavg..F. - 32)
######################################################### END CLIMATE INPUTS ####################################################################
######################################################### CALCULATE WB VARIABLES ################################################################
AllDailyWB<-list()
DailyWB<-Historical
for(i in 1:nrow(wb_sites)){
ID = wb_sites$SiteID[i] # KW: was previously "Site"
Lat = wb_sites$Lat[i]
Lon = wb_sites$Lon[i]
Elev = wb_sites$Elev[i] # KW: was previously "Elevation"
Aspect = wb_sites$Aspect[i]
Slope = wb_sites$Slope[i]
SWC.Max = wb_sites$SWC.Max[i]
Wind = wb_sites$Wind[i]
Snowpack.Init = wb_sites$Snowpack.Init[i]
Soil.Init = wb_sites$Soil.Init[i]
Shade.Coeff = wb_sites$Shade.Coeff[i]
#Calculate daily water balance variables
DailyWB$ID = ID
DailyWB$doy <- yday(DailyWB$Date)
DailyWB$daylength = get_daylength(DailyWB$Date, Lat)
DailyWB$jtemp = as.numeric(get_jtemp(Lon, Lat))
DailyWB$F = get_freeze(DailyWB$jtemp, DailyWB$tmean_C)
DailyWB$RAIN = get_rain(DailyWB$ppt_mm, DailyWB$F)
DailyWB$SNOW = get_snow(DailyWB$ppt_mm, DailyWB$F)
DailyWB$MELT = get_melt(DailyWB$tmean_C, DailyWB$jtemp, hock=4, DailyWB$SNOW, Snowpack.Init)
DailyWB$PACK = get_snowpack(DailyWB$jtemp, DailyWB$SNOW, DailyWB$MELT)
DailyWB$W = DailyWB$MELT + DailyWB$RAIN
if(Method == "Hamon"){
DailyWB$PET = ET_Hamon_daily(DailyWB)
} else {
if(Method == "Penman-Monteith"){
DailyWB$PET = ET_PenmanMonteith_daily(DailyWB)
} else {
if(Method == "Oudin"){
DailyWB$PET = get_OudinPET(DailyWB$doy, Lat, DailyWB$PACK, DailyWB$tmean_C, Slope, Aspect, Shade.Coeff)
} else {
print("Error - PET method not found")
}
}
}
DailyWB$PET = modify_PET(DailyWB$PET, Slope, Aspect, Lat, Shade.Coeff)
DailyWB$W_PET = DailyWB$W - DailyWB$PET
DailyWB$SOIL = get_soil(DailyWB$W, Soil.Init, DailyWB$PET, DailyWB$W_PET, SWC.Max)
DailyWB$DSOIL = diff(c(Soil.Init, DailyWB$SOIL))
DailyWB$AET = get_AET(DailyWB$W, DailyWB$PET, DailyWB$SOIL, Soil.Init)
DailyWB$W_ET_DSOIL = DailyWB$W - DailyWB$AET - DailyWB$DSOIL
DailyWB$D = DailyWB$PET - DailyWB$AET
DailyWB$GDD = get_GDD(DailyWB$tmean_C, T.Base)
AllDailyWB[[i]] = DailyWB
}
WBData<-do.call(rbind,AllDailyWB)
Lastly, we aggregate the 10 WBM outputs by date and grab the means for a handful of the WBM variables we want to compare. We also convert the output (mm) to inches.
WBData_coded <- do.call(rbind, AllDailyWB) %>%
group_by(Date, GCM) %>%
# convert to inches
dplyr::summarize(mean_runoff = mean(W_ET_DSOIL, na.rm = TRUE)/25.4,
mean_pet = mean(PET, na.rm = TRUE)/25.4,
mean_rain = mean(RAIN, na.rm = TRUE)/25.4,
mean_snow = mean(SNOW, na.rm = TRUE)/25.4,
mean_melt = mean(MELT, na.rm = TRUE)/25.4,
mean_snowpack = mean(PACK, na.rm = TRUE)/25.4,
mean_aet = mean(AET, na.rm = TRUE)/25.4) %>%
# Remove "warm up" year
filter(year(as_date(Date)) > 1979)
Next, we can pull in and organize the data associated with the other data products on AWS: the gridded product, and the centroid product.
Here is the centroid data downloaded directly from AWS for PETE (https://parkfutures.s3.us-west-2.amazonaws.com/park-centroid-wb/index.html). We also remove the first year of data (1979), since we are under the impression that this year is used as the “warm up” period, and because our gridded data set only has data starting in 1980. Data are reported in inches in this data set. We are also under the impression that the Oudin method was selected for calculating PET.
WBData_online_centroid <- read_csv("aws_data_downloads/centroid_data/PETE_water_balance_historical.csv") %>%
dplyr::filter(year(as_date(Date)) > 1979)
Here, we pull in the gridded data, which was downloaded directly from AWS, then clipped to PETE’s park boundary. Gridded data starts in 1980 (which could be explained by the fact that 1979 is used as a “warm up”). We are also under the impression that the Oudin method was selected for calculating PET. The python script provided to us that represents the workflow used to develop the gridded dataset can be found at “aws_python_script/tercek_script.py”.
WBData_online_grid <- rast("aws_data_downloads/gridded_data/PETE_runoff_historic_gridmet_1980_2023.tif")
# next we clip the grid to only grids that intersect our wb_sites points in PETE:
NPS_points <- st_transform(sf::st_as_sf(wb_sites, coords = c("Lon", "Lat"), crs = 4326), st_crs(WBData_online_grid))
grid_runoff <- terra::extract(WBData_online_grid, NPS_points) %>%
pivot_longer(cols = -ID) %>%
dplyr::group_by(name) %>%
# convert to inches
dplyr::summarize(mean_runoff = mean(value, na.rm = TRUE)/25.4/10) %>%
filter(name != "runoff_NA") %>%
mutate(date = as_date(sub(".*_", "", name))) %>%
filter(year(date) < 2023)
WBData_online_grid <- rast("aws_data_downloads/gridded_data/PETE_pet_historic_gridmet_1980_2023.tif")
grid_pet <- terra::extract(WBData_online_grid, NPS_points) %>%
pivot_longer(cols = -ID) %>%
dplyr::group_by(name) %>%
# convert to inches
dplyr::summarize(mean_pet = mean(value, na.rm = TRUE)/25.4/10) %>%
filter(name != "PET_NA") %>%
mutate(date = as_date(sub(".*_", "", name))) %>%
filter(year(date) < 2023)
WBData_online_grid <- rast("aws_data_downloads/gridded_data/PETE_rain_historic_gridmet_1980_2023.tif")
grid_rain <- terra::extract(WBData_online_grid, NPS_points) %>%
pivot_longer(cols = -ID) %>%
dplyr::group_by(name) %>%
# convert to inches
dplyr::summarize(mean_rain = mean(value, na.rm = TRUE)/25.4/10) %>%
filter(name != "rain_NA") %>%
mutate(date = as_date(sub(".*_", "", name))) %>%
filter(year(date) < 2023)
Here I am plotting each of the different WBM runoff products against each other for comparison.
year_colors <- c("#002EA3", "#C3F73A", "#E70870", "#256BF5", "#745CFB", "#FFCA3A", "#578010", "#56104E")
interpolated_colors <- colorRampPalette(year_colors)(length(unique(year(WBData_coded$Date))))
one <- ggplot() +
ggtitle("Runoff in") +
theme_bw() +
xlab("Generated Data from Code") +
ylab("Centroid Data on AWS") +
geom_point(aes(x = WBData_coded$mean_runoff,
y = WBData_online_centroid$`runoff in`,
color = factor(year(WBData_coded$Date))),
size = 1) +
scale_color_manual(values = interpolated_colors) +
labs(color = "Year")
two <- ggplot() +
theme_bw() +
xlab("Generated Data from Code") +
ylab("Gridded Data on AWS") +
geom_point(aes(x = WBData_coded$mean_runoff,
y = grid_runoff$mean_runoff,
color = factor(year(WBData_coded$Date))),
size = 1) +
scale_color_manual(values = interpolated_colors) +
labs(color = "Year")
three <- ggplot() +
theme_bw() +
xlab("Centroid Data on AWS") +
ylab("Gridded Data on AWS") +
geom_point(aes(x = WBData_online_centroid$`runoff in`,
y = grid_runoff$mean_runoff,
color = factor(year(WBData_coded$Date))),
size = 1) +
scale_color_manual(values = interpolated_colors) +
labs(color = "Year")
ggpubr::ggarrange(one + theme(legend.position = "none"),
two + theme(legend.position = "none"),
three + theme(legend.position = "none"),
ncol = 1, nrow = 3,
common.legend = TRUE, legend = "right")
All three runoff products plotted together through time:
ggplotly(
ggplot() +
geom_point(data = WBData_online_centroid, aes(x = as_date(Date), y = `runoff in`), color = "#E70870", alpha = 1) +
geom_point(data = grid_runoff, aes(x = date, y = mean_runoff), color = "#256BF5", alpha = 0.6) +
geom_line(data = WBData_coded, aes(x = as_date(Date), y = mean_runoff), color = "black", cex = 0.2) +
theme_bw() +
labs(x = "Date", y = "Runoff in")
# geom_text(data = WBData_online_centroid[1, ], aes(x = as.Date(Date), y = `runoff in`), label = "Centroid Data on AWS", color = "#E70870", vjust = -35, hjust = -1) +
# geom_text(data = grid_runoff[1, ], aes(x = date, y = mean_runoff), label = "Gridded Data on AWS", color = "#256BF5", vjust = -33, hjust = -1) +
# geom_text(data = WBData_coded[1, ], aes(x = as.Date(Date), y = mean_runoff), label = "Generated Data from Code", color = "black", vjust = -31, hjust = -0.8)
)
Figure Legend: Centroid Data on AWS, Generated Data from Code, Gridded Data o AWS.
Here I am plotting each of the different WBM PET products against each other for comparison.
one <- ggplot() +
ggtitle("PET in") +
theme_bw() +
xlab("Generated Data from Code") +
ylab("Centroid Data on AWS") +
geom_point(aes(x = WBData_coded$mean_pet,
y = WBData_online_centroid$`PET in`,
color = factor(year(WBData_coded$Date))),
size = 1) +
scale_color_manual(values = interpolated_colors) +
labs(color = "Year")
two <- ggplot() +
theme_bw() +
xlab("Generated Data from Code") +
ylab("Gridded Data on AWS") +
geom_point(aes(x = WBData_coded$mean_pet,
y = grid_pet$mean_pet,
color = factor(year(WBData_coded$Date))),
size = 1) +
scale_color_manual(values = interpolated_colors) +
labs(color = "Year")
three <- ggplot() +
theme_bw() +
xlab("Centroid Data on AWS") +
ylab("Gridded Data on AWS") +
geom_point(aes(x = WBData_online_centroid$`PET in`,
y = grid_pet$mean_pet,
color = factor(year(WBData_coded$Date))),
size = 1) +
scale_color_manual(values = interpolated_colors) +
labs(color = "Year")
ggpubr::ggarrange(one + theme(legend.position = "none"),
two + theme(legend.position = "none"),
three + theme(legend.position = "none"),
ncol = 1, nrow = 3,
common.legend = TRUE, legend = "right")
All three PET products plotted together through time:
ggplotly(
ggplot() +
geom_point(data = WBData_online_centroid, aes(x = as_date(Date), y = `PET in`), color = "#E70870", alpha = 1) +
geom_point(data = grid_pet, aes(x = date, y = mean_pet), color = "#256BF5", alpha = 0.6) +
geom_line(data = WBData_coded, aes(x = as_date(Date), y = mean_pet), color = "black", cex = 0.2) +
theme_bw() +
labs(x = "Date", y = "PET in")
)
Figure Legend: Centroid Data on AWS, Generated Data from Code, Gridded Data on AWS.
Here I am plotting each of the different WBM rain products against each other for comparison.
one <- ggplot() +
ggtitle("Rain in") +
theme_bw() +
xlab("Generated Data from Code") +
ylab("Centroid Data on AWS") +
geom_point(aes(x = WBData_coded$mean_rain,
y = WBData_online_centroid$`rain in`,
color = factor(year(WBData_coded$Date))),
size = 1) +
scale_color_manual(values = interpolated_colors) +
labs(color = "Year")
two <- ggplot() +
theme_bw() +
xlab("Generated Data from Code") +
ylab("Gridded Data on AWS") +
geom_point(aes(x = WBData_coded$mean_rain,
y = grid_rain$mean_rain,
color = factor(year(WBData_coded$Date))),
size = 1) +
scale_color_manual(values = interpolated_colors) +
labs(color = "Year")
three <- ggplot() +
theme_bw() +
xlab("Centroid Data on AWS") +
ylab("Gridded Data on AWS") +
geom_point(aes(x = WBData_online_centroid$`rain in`,
y = grid_rain$mean_rain,
color = factor(year(WBData_coded$Date))),
size = 1) +
scale_color_manual(values = interpolated_colors) +
labs(color = "Year")
ggpubr::ggarrange(one + theme(legend.position = "none"),
two + theme(legend.position = "none"),
three + theme(legend.position = "none"),
ncol = 1, nrow = 3,
common.legend = TRUE, legend = "right")
All three rain products plotted together through time:
ggplotly(
ggplot() +
geom_point(data = WBData_online_centroid, aes(x = as_date(Date), y = `rain in`), color = "#E70870", alpha = 1) +
geom_point(data = grid_rain, aes(x = date, y = mean_rain), color = "#256BF5", alpha = 0.6) +
geom_line(data = WBData_coded, aes(x = as_date(Date), y = mean_rain), color = "black", cex = 0.2) +
theme_bw() +
labs(x = "Date", y = "Rain in")
)
Figure Legend: Centroid Data on AWS, Generated Data from Code, Gridded Data o AWS.
Based on discussions with the NPS group, we understand that the published centroid data and the gridded data products were developed with different intentions, and subsequently, are different by design. For example:
The gridded dataset’s workfow has the soil water storage start at full capacity, whereas the WBM R code has the soil water storage start at 0.
The WBM R code’s calculation for Oudin PET is slightly different: it incorporates a shade coefficient to reduce PET (though it is set to NULL and therefore not in use), whereas the gridded workflow does not include this option. The WBM R code also sets PET to 0 if there is at least 2 mm of snowpack, whereas the gridded product sets PET to 0 if there is any snow (snowpack > 0).
However, we are under the impression that the park centroid data was generated using the same workflow outlined in the WBM R code. Though they are incredibly similar, they are not identical. Here, we lay out a few potential explanations for the differences in their outputs.
get_OudinPET() incorporates aspects of
modify_PET(). Both are used in the published tutorial on
AWS in generating WBM variables. (See lines 107-120 above.) Below are
the two functions:get_OudinPET = function(doy, lat, snowpack, tmean, slope, aspect, shade.coeff=NULL){
d.r = 1 + 0.033*cos((2*pi/365)*doy)
declin = 0.409*sin((((2*pi)/365)*doy)-1.39)
lat.rad = (pi/180)*lat
sunset.ang = acos(-tan(lat.rad)*tan(declin))
R.a = ((24*60)/pi)*0.082*d.r*((sunset.ang*sin(lat.rad)*sin(declin)) + (cos(lat.rad)*cos(declin)*sin(sunset.ang)))
Oudin = ifelse(snowpack>2,0,ifelse(tmean>-5,(R.a*(tmean+5)*0.408)/100,0))
Folded_aspect = abs(180-abs((aspect)-225))
Heatload = (0.339+0.808*cos(lat*(pi/180))*cos(slope*(pi/180)))-(0.196*sin(lat.rad)*sin(slope*(pi/180)))-(0.482*cos(Folded_aspect*(pi/180))*sin(slope*(pi/180)))
sc = ifelse(!is.null(shade.coeff), shade.coeff, 1)
OudinPET = Oudin * Heatload * sc
return(OudinPET)
}
modify_PET = function(pet, slope, aspect, lat, freeze, shade.coeff=NULL){
f.aspect = abs(180 - abs(aspect - 225))
lat.rad = ifelse(lat > 66.7, (66.7/180)*pi, (lat/180)*pi)
slope.rad = (slope/180)*pi
aspect.rad = (f.aspect/180)*pi
heat.load = 0.339+0.808*cos(lat.rad)*cos(slope.rad) - 0.196*sin(lat.rad)*sin(slope.rad) - 0.482*cos(aspect.rad)*sin(slope.rad)
sc = ifelse(!is.null(shade.coeff), shade.coeff, 1)
freeze = ifelse(freeze == 0,0,1)
PET.Lutz = pet*heat.load*sc*freeze
return(PET.Lutz)
}
Of note, the exclusion of the modify_PET() step does not
lead to a major change in values. Below is a plot of the difference
between a runoff output that uses modify_PET() and the
runoff output that does not use modify_PET():
n <- nrow(wb_sites)
#Threshold temperature (deg C) for growing degree-days calculation
T.Base = 0
#Method for PET calculation
Method = "Oudin" #Oudin is default method for daily PRISM and MACA data (containing only Tmax, Tmin, and Date).
#Date format
DateFormat = "%m/%d/%Y"
############################################################ CREATE CLIMATE INPUTS #############################################################
#### Historical
# Convert pr.In to mm and F to C
Historical$ppt_mm <- (Historical$Precip..in.*25.4)
Historical$tmax_C <- 5/9*(Historical$Tmax..F. - 32)
Historical$tmin_C <- 5/9*(Historical$Tmin..F. - 32)
Historical$tmean_C <- 5/9*(Historical$Tavg..F. - 32)
######################################################### END CLIMATE INPUTS ####################################################################
######################################################### CALCULATE WB VARIABLES ################################################################
AllDailyWB<-list()
DailyWB<-Historical
for(i in 1:nrow(wb_sites)){
ID = wb_sites$SiteID[i] # KW: was previously "Site"
Lat = wb_sites$Lat[i]
Lon = wb_sites$Lon[i]
Elev = wb_sites$Elev[i] # KW: was previously "Elevation"
Aspect = wb_sites$Aspect[i]
Slope = wb_sites$Slope[i]
SWC.Max = wb_sites$SWC.Max[i]
Wind = wb_sites$Wind[i]
Snowpack.Init = wb_sites$Snowpack.Init[i]
Soil.Init = wb_sites$Soil.Init[i]
Shade.Coeff = wb_sites$Shade.Coeff[i]
#Calculate daily water balance variables
DailyWB$ID = ID
DailyWB$doy <- yday(DailyWB$Date)
DailyWB$daylength = get_daylength(DailyWB$Date, Lat)
DailyWB$jtemp = as.numeric(get_jtemp(Lon, Lat))
DailyWB$F = get_freeze(DailyWB$jtemp, DailyWB$tmean_C)
DailyWB$RAIN = get_rain(DailyWB$ppt_mm, DailyWB$F)
DailyWB$SNOW = get_snow(DailyWB$ppt_mm, DailyWB$F)
DailyWB$MELT = get_melt(DailyWB$tmean_C, DailyWB$jtemp, hock=4, DailyWB$SNOW, Snowpack.Init)
DailyWB$PACK = get_snowpack(DailyWB$jtemp, DailyWB$SNOW, DailyWB$MELT)
DailyWB$W = DailyWB$MELT + DailyWB$RAIN
if(Method == "Hamon"){
DailyWB$PET = ET_Hamon_daily(DailyWB)
} else {
if(Method == "Penman-Monteith"){
DailyWB$PET = ET_PenmanMonteith_daily(DailyWB)
} else {
if(Method == "Oudin"){
DailyWB$PET = get_OudinPET(DailyWB$doy, Lat, DailyWB$PACK, DailyWB$tmean_C, Slope, Aspect, Shade.Coeff)
} else {
print("Error - PET method not found")
}
}
}
#DailyWB$PET = modify_PET(DailyWB$PET, Slope, Aspect, Lat, Shade.Coeff)
DailyWB$W_PET = DailyWB$W - DailyWB$PET
DailyWB$SOIL = get_soil(DailyWB$W, Soil.Init, DailyWB$PET, DailyWB$W_PET, SWC.Max)
DailyWB$DSOIL = diff(c(Soil.Init, DailyWB$SOIL))
DailyWB$AET = get_AET(DailyWB$W, DailyWB$PET, DailyWB$SOIL, Soil.Init)
DailyWB$W_ET_DSOIL = DailyWB$W - DailyWB$AET - DailyWB$DSOIL
DailyWB$D = DailyWB$PET - DailyWB$AET
DailyWB$GDD = get_GDD(DailyWB$tmean_C, T.Base)
AllDailyWB[[i]] = DailyWB
}
WBData_coded_no_modify_PET <- do.call(rbind, AllDailyWB) %>%
group_by(Date, GCM) %>%
# convert to inches
dplyr::summarize(mean_runoff = mean(W_ET_DSOIL, na.rm = TRUE)/25.4,
mean_pet = mean(PET, na.rm = TRUE)/25.4,
mean_rain = mean(RAIN, na.rm = TRUE)/25.4,
mean_snow = mean(SNOW, na.rm = TRUE)/25.4,
mean_melt = mean(MELT, na.rm = TRUE)/25.4,
mean_snowpack = mean(PACK, na.rm = TRUE)/25.4,
mean_aet = mean(AET, na.rm = TRUE)/25.4) %>%
# Remove "warm up" year
filter(year(as_date(Date)) > 1979)
ggplot() +
#geom_point(data = WBData_coded, aes(x = as_date(Date), y = mean_runoff), color = "black", cex = 1.5) +
geom_col(aes(x = as_date(WBData_coded$Date), y = WBData_coded$mean_runoff - WBData_coded_no_modify_PET$mean_runoff), color = "#578010") +
theme_bw() +
ggtitle("Difference in Runoff (inches)") +
xlab("Date") +
ylab("With - Without `modify_PET()`")
Though we assumed that the coordinates provided in the folder hosted on AWS (“aws_r_tutorial/PETE_site_parameters.csv”) were used in developing PETE’s published centroid data, it is possible that this is not the case.
We were unable to find the exact historical climate dataset used in the AWS tutorial. Instead, we used PETE’s centroid’s historical GridMET data posted on AWS (https://parkfutures.s3.us-west-2.amazonaws.com/maca-tprh-data/index.html). If a different data set was used to generate the centroid data, this could also explain the discrepancies.
The code we are running starts the WBM in the year 1979. Though we believe that 1979 is also the year used to generate the centroid data, it is possible another year was used as the starting point.